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[i]  A  new  code  for  solving  radiation  belt  diffusion  equations  has  been  developed  and 
applied  to  the  2-D  bounce-averaged  energy  pitch  angle  quasi-linear  diffusion  equation. 

The  code  uses  Monte  Carlo  methods  to  solve  Ito  stochastic  differential  equations 
(SDEs)  which  are  mathematically  equivalent  to  radiation  belt  diffusion  equations.  We 
show  that  our  SDE  code  solves  the  diffusion  equation  with  off-diagonal  diffusion 
coefficients  in  contrast  to  standard  finite  difference  codes  which  are  generally  unstable 
when  off-diagonal  diffusion  coefficients  are  included.  Our  results  are  in  excellent 
agreement  with  previous  results.  We  have  also  investigated  effects  of  assuming  purely 
parallel  propagating  electromagnetic  waves  when  calculating  the  diffusion  coefficients  and 
find  that  this  assumption  leads  to  errors  of  more  than  an  order  of  magnitude  in  flux  at 
some  equatorial  pitch  angles  for  the  specific  chorus  wave  model  we  use.  Further  work  is 
needed  to  investigate  the  sensitivity  of  our  results  to  the  wave  model  parameters. 

Generalization  of  the  method  to  3-D  is  straightforward,  thus  making  this  method  a  very 
promising  new  way  to  investigate  the  relative  roles  of  pitch  angle,  energy,  and  radial 
diffusion  in  radiation  belt  dynamics. 

Citation:  Tao,  X.,  A.  A.  Chan,  J.  M.  Albert,  and  J.  A.  Miller  (2008),  Slochaslic  modeling  of  mullidimcnsional  diffusion  in  the 
radiation  belts,  J.  Geophys.  Res.,  113,  A07212,  doi:10.1029/2007JA0 12985. 


1.  Introduction 

[2]  The  Earth’s  outer  radiation  belt  is  very  dynamic,  and 
electron  fluxes  can  vary  by  several  orders  of  magnitude 
during  storm  times,  which  makes  it  very  hazardous 
to  spacecrafts  and  astronauts  [e.g.,  Baker  et  ai,  1997]. 
Quasi-linear  diffusion  theory  has  been  used  to  evaluate 
dynamic  changes  of  particle  fluxes  in  the  radiation  belts 
[Albert,  2004;  Albert  and  Young,  2005;  Horne  and  Thorne, 
2003;  Horne  et  ai,  2003].  Using  the  quasi-linear  diffusion 
theory  to  model  radiation  belt  dynamics  requires  at  least  two 
kinds  of  computations:  numerical  solution  of  a  diffusion 
equation,  which  is  a  one-dimensional  or  multidimensional 
Fokkcr-Planck  equation,  depending  on  diffusion  processes 
we  are  interested  in,  and  calculation  of  diffusion  coefficients. 

[3]  Albert  [2004]  has  shown  that  numerical  problems 
arise  when  applying  standard  finite  difference  methods  to 
pitch  angle  and  energy  diffusion  equations  because  of 
rapidly  varying  off-diagonal  diffusion  coefficients.  Albert 
and  Young  [2005]  developed  a  method  for  the  2-D  diffusion 
equation  which  diagonalizes  the  diffusion  tensor  by  trans¬ 
forming  to  a  new  set  of  coordinates  and  solves  the  trans¬ 
formed  equation  by  simple  finite  difference  methods.  In  this 
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woric  we  introduce  another  method  which  uses  probabilistic 
representations  of  solutions  of  Fokkcr-Planck  equations 
[Freidlin,  1985;  Costantini  et  ai,  1998]  via  stochastic 
differential  equations  (SDEs),  and  we  develop  a  2-D  code 
for  solving  pitch  angle  and  energy  diffusion  equations. 
Compared  with  finite  difference  methods,  the  SDE  method 
has  three  main  advantages.  First,  the  SDE  method  is  very 
efficient  when  solutions  on  only  a  small  number  of  points 
are  desired,  particularly  when  applied  to  high-dimensional 
problems,  and  it  is  easy  to  code  and  parallelize,  with 
parallelization  efficiency  close  to  one.  Second,  with  the 
SDE  method,  wc  are  able  to  handle  complicated  boundary 
geometry  other  than  constant-coordinate  boundaries  (see 
section  2.2).  Third,  generalization  of  the  SDE  method  to 
higher  dimensions  is  straightforward,  and  we  expect  the 
method  to  be  applicable  to  general  3-D  radiation  belt 
diffusion  equations.  For  more  applications  of  similar 
methods  using  relations  between  Fokkcr-Planck  equations 
and  SDEs,  see,  e.g.,  Zhang  [1999],  Albright  et  ai  [2003], 
Alanko-Huotari  et  ai  [2007],  Qin  et  ai  [2005],  and  Yamada 
etai  [1998]. 

[4]  Besides  solving  diffusion  equations,  correctly  calcu¬ 
lating  quasi-linear  diffusion  coefficients  is  also  important 
for  numerical  modeling  of  the  radiation  belt  dynamics  using 
quasi-linear  theory.  Albert  [2005]  and  Glauert  and  Horne 
[2005]  have  shown  full  calculations  of  diffusion  coefficients 
for  cyclotron  resonant  wave-particle  interactions,  where  up 
to  ^  =  ±  5  resonances  are  included.  However,  the  full 
calculation  of  diffusion  coefficients  is  very  time  consuming. 
Summers  [2005]  derived  simplified  formulae  for  coeffi¬ 
cients  with  a  parallel  propagation  approximation  (and  hence 
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only  the  «  =  -1  resonance  is  included  [Albert,  2007]),  and 
the  computation  becomes  much  fasten  Shprits  et  al  [2006] 
calculated  bounce-averaged  pitch  angle  and  energy  diffu¬ 
sion  eoeffieients  D^^oao  ^pp  with  the  parallel  propaga¬ 
tion  approximation  for  £  <  1  MeV  partieles  and  compared 
them  with  fully  calculated  coefficients  from  the  PADIE  code 
of  Glauert  and  Horne  [2005].  They  concluded  that  coef¬ 
ficients  for  field-aligned  waves  are  elose  to  coefficients  for 
waves  with  mildly  oblique  wave  normal  angle  distribution 
from  the  PADIE  code.  However,  using  the  wave  model 
from  Horne  et  al.  [2005],  we  compute  particle  fluxes  and 
wc  show  that  for  £"  =  2  MeV  electrons,  and  Dpp 

calculated  with  the  parallel  propagation  approximation 
produce  flux  differences  of  about  1  order  of  magnitude  at 
some  pitch  angles,  compared  to  using  fully  calculated 
coefficients.  Furthermore,  we  show  that  by  including  off- 
diagonal  terms  in  the  calculation,  the  parallel  propagation 
approximation  also  produces  large  errors  in  fluxes  for  both 
£■  =  0.5  MeV  and  2  MeV  electrons  at  small  pitch  angles. 

[5]  The  remainder  of  this  paper  is  organized  as  follows. 
The  SDE  method  and  its  numerical  implementation  are 
introduced  in  section  2.  In  section  3  we  present  the 
application  of  the  SDE  method  to  a  bounce-averaged 
radiation  belt  pitch  angle  and  energy  diffusion  equation. 
After  describing  the  implementation  of  the  SDE  method  for 
the  pitch  angle  energy  equation  (section  3.1),  we  show 
comparisons  between  results  from  the  SDE  method  and  the 
Albert  and  Young  [2005]  transformation  method  (section  3.2). 
Then  fluxes  calculated  from  diffusion  coefficients  with  the 
parallel  propagation  approximation  [Summers,  2005]  are 
compared  with  fluxes  computed  with  coefficients  from  full 
quasi-linear  theory  [Albert,  2005]  (section  3.3).  We  summarize 
our  work  and  discuss  future  work  in  section  4. 

2.  SDE  Method 


Gaussian  random  number.  The  m  vector  b  and  the  m  x  n 
matrix  a  are  coefficients  that  determine  the  values  of  X(0,  they 
will  be  directly  related  to  the  coefficients  of  a  corresponding 
diffusion  equation  in  section  2.2.  Stepping  equation  (1 )  in  time 
generates  a  random  walk  trajectory  through  X  space. 

[8]  Note  that  SDEs  may  be  formulated  using  two  main 
mathematical  methods:  the  Ito  method  and  the  Stratonovich 
method  [Gardiner,  1985].  In  this  work  we  use  Ito  SDEs 
because  they  are  directly  related  to  diffusion  equations  of 
interest  for  the  radiation  belts  and  they  are  mathematically 
more  convenient  [Oksendal,  \992\  Freidlin,  1 985;  Costantini 
et  al.,  1998]. 

2.2.  Probabilistic  Representation  of  Solutions  of 
Diffusion  Equations 

[9]  To  solve  a  diffusion  equation  using  SDEs,  wc  can  first 
write  the  diffusion  equation  in  Fokker-Planck  form  and  then 
obtain  equivalent  “time-forward”  SDEs  from  the  diffusion 
equation.  These  time-forward  SDEs  can  then  be  used  to 
simulate  particle  trajectories  using  a  Monte  Carlo  technique, 
and  the  distribution  of  particles  at  any  given  time  ean  be 
obtained  by  binning  partieles  in  phase  space.  This  time- 
forward  SDE  method  is  presented  in  Appendix  A  to  show 
local  effects  of  off-diagonal  terms  on  the  distribution  of 
particles.  Alternatively,  in  this  section  we  present  a  “time- 
backward”  SDE  method,  where  solutions  of  diffusion 
equations  are  represented  by  the  mean  value  of  a  functional 
of  trajectories  of  a  stochastic  process  [Freidlin,  1985].  This 
is  the  method  used  in  our  current  SDE  code.  Compared  with 
the  time-forward  method,  the  time-backward  method  is 
more  efficient  when  solutions  on  fewer  points  are  of 
interest,  and  it  is  better  for  handling  a  variety  of  boundary 
conditions. 

[10]  To  introduce  the  time-baekward  SDE  method,  let  us 
first  consider  a  ^/-dimensional  diffusion  equation  written  as 


[6]  Our  SDE  code  is  based  on  mathematical  results  which 
show  that  solutions  of  diffusion  equations  can  be  obtained 
using  an  equivalent  stochastic  process.  Thus,  we  first  give  a 
description  of  a  stochastic  process  using  Ito  stochastic 
differential  equations  in  section  2.1.  Then  we  show  how 
these  lead  to  probabilistic  representations  of  solutions  of 
diffusion  equations  in  section  2.2. 

2.1.  Ito  Stochastic  DifTerential  Equations 

[7]  Stochastic  differential  equations  (SDEs)  are  used  to 
describe  stochastic  processes.  They  differ  from  ordinary 
differential  equations  by  having  terms  involving  random 
variables  [Gardiner,  1985;  Freidlin,  1985].  A  general 
w-dimensional  SDE  with  an  /i-dimensional  Wiener  process 
is  written  as 


dXit)  =  b(X.  t)  dt  +  <t(X,  t)  ^/ W(t),  ( 1 ) 

where  the  m  vector  X  represents  an  w-dimensional 
stochastic  process  {X\,  Y2,. . .,  Y^).  Throughout  this  work, 
stochastic  processes  are  indicated  by  uppercase  characters, 
and  their  values  at  a  given  time  are  represented  by 
corresponding  lowercase  characters.  The  n  vector  W  is 
an  «-dimensional  Wiener  process  {W^,  W2,  . . W„)  and 
i/W(0  =  W(r  +  dt)  -  W(0  [Gardiner,  1985];  an  increment  of 
a  one-dimensional  Wiener  process  is  proportional  to  a 


Of 


dt  2 

d 

E 

/=i 


OxiOxj 

Of 


{t.x) 


+  E  of. 


with  initial  and  Dirichlet  boundary  conditions 
/(0,x)  =go(x),  xeD. 


(2) 

(3) 


/(Lx)  =gi(/.x),  xeOD.  (4) 

Here  D  is  the  domain  of  the  problem  with  boundary  OD,  and 
gi(0,  x)  =  go(x)  on  OD.  Note  that  OD  is  not  restricted  to 
constant  coordinate  surfaces  in  the  SDE  method  [Freidlin, 
1985]. 

[11]  The  solution  /(x,  t)  of  equation  (2)  is  related  to  the 
following  ^/-dimensional  stochastic  process: 

dX(s)  =  b(t  -  s,X)  ds  ^  (T(t  ~  s,X)  dW(s),  0  <  s  <  t,  (5) 

where  X(s  =  0)  =  x  and  W(.s)  is  a  ^/-dimensional  Wiener 
process.  Here  the  d  x  d  matrix  a  is  defined  by  crer^  =  a. 
Note  that  a  is  not  uniquely  determined  by  this  equation,  but 
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according  to  Levy’s  theorem  [Zhang,  1999;  Freidlin,  1985], 
different  choiees  of  a  generate  equivalent  stochastic 
processes  that  yield  the  same  solution  of  the  diffusion 
equation  (2).  Also,  note  that  equation  (5)  is  a  time-backward 
SDE:  at  5  =  0,  we  evaluate  b  and  cr  at  time  t,  while  at  5  - 
we  evaluate  b  and  a  at  time  zero.  The  solution /(x,  /)  is  then 
represented  by  the  stochastic  process  defined  in  equation  (5) 
as 


pared  with  results  from  Albert  and  Young  [2005]  to  show  that 
the  SDE  code  is  eapable  of  solving  the  diffusion  equation 
with  off-diagonal  diffusion  coefficients.  To  show  the  effect  of 
diffusion  coefficients  with  the  parallel  propagation  approxi¬ 
mation  [Summers,  2005]  on  particle  fluxes,  we  solve  the 
diffusion  equation  using  these  diffusion  coefficients  and  in 
section  3.3  the  results  arc  compared  with  those  obtained  from 
fiilly  calculated  coefficients. 


/{X,/)  =  E{F,),  (6) 

where  E  denotes  the  expectation  value  and  F,  is  defined  by 

F.  =  {  (7) 

(  gi  (;  -  r,X|,,^)  exp(y'|,^^),  r  <  ;, 

where  r  has  the  value  of  s  when  the  stochastic  process  Xf?) 
exits  from  the  boundary  OD  for  the  first  time  and  T(5)  is 
defined  by 


Y{s)  =  r  c{t-r.X{r))dr.  (8) 

Jo 

[12]  Numerical  calculation  of/ can  be  constructed  easily 
from  equations  (6)  (8).  To  obtain  /(x,  /),  we  sample  a 
number  of  trajectories  of  the  stochastic  process  defined  by 
equation  (5)  starting  from  x  and  s  =  0,  using  a  Monte  Carlo 
technique.  The  simulation  of  a  trajectory  will  stop  either  by 
reaching  the  initial  condition  at  s  =  t  (where  time  =  0)  or  by 
reaching  the  boundary  of  the  domain  D  at  s  =  r,  whichever 
comes  first,  and  returns  a  value  defined  by  equation  (7). 
Then  we  use  the  average  of  values  returned  by  all  trajecto¬ 
ries  to  approximate  /(x,  /).  This  process  is  repeated  if  we 
want  to  calculate  /at  other  points. 

[13]  Now  let  us  also  consider  a  particular  type  of 
Neumann  boundary  condition  that  is  commonly  encoun¬ 
tered  in  radiation  belt  diffusion  equations: 

V/.n  =  0.  x€<9,D,  (9) 

where  V/=  {df/0x\  df/Ox^,. . .,  df/0)^),  the  boundary  OyD 
is  the  part  of  OD  with  the  Neumann  condition,  and  n  is  the 
inward  unit  normal  vector  on  d\D.  General  methods  for 
implementing  Neumann  boundary  conditions  in  SDE  solu¬ 
tions  are  given  by  Freidlin  [1985]  and  Costantini  et  al. 
[1998];  here  we  simply  note  that  condition  (9)  ean  be 
enforced  in  our  numerical  calculation  of /(x,  /)  as  follows: 
Every  time  a  trajectory  reaches  the  Neumann  boundary  0\D, 
we  immediately  reflect  it  about  the  normal  vector  n  et 

al.,  2004].  This  trajectory  will  later  be  stopped  by  either 
reaching  the  initial  condition  or  a  Dirichlet  boundary,  and  at 
that  time  the  trajectory  returns  a  value  defined  by  equation  (7). 

3.  Application 

[14]  In  this  section,  we  apply  the  above  (see  section  2)  SDE 
method  to  a  bounce-averaged  pitch  angle  and  energy  diffu¬ 
sion  equation  [Albert,  2004].  In  section  3.1  we  derive  the 
stochastic  process  used  to  solve  the  diffusion  equation.  In 
section  3.2  fluxes  calculated  using  the  SDE  code  are  com- 


3.1.  Application  to  Pitch  Angle  and  Energy 
Diffusion  Equations 

[15]  We  apply  the  above  SDE  method  to  the  boiinec- 
averaged  pitch  angle  and  energy  diffusion  equation  written 
in  equatorial  pitch  angle  and  momentum  (oo,  p) 


(If  1  (1  1  Of  ^  (lf\ 

0l  ~  Gp  dna  ^  P  Oon  Op) 

\  0  (  \  Of  0f\ 

G  ^  pO^o"  ^  j  ' 


(10) 


where  and  D^p  arc  bounce-averaged  pitch 

angle,  mixed,  and  momentum  diffusion  eoeffieients  [Albert, 
2004],  Here  G  is  a  Jacobian  factor,  G=p^7\ao)  sin(ao)  cosftvo), 
and  r(ao)  ^  1.30  -  0.56  sin(ao)  is  the  normalized  bounce 
period.  Initial  and  boundary  conditions  are  chosen  to  be  the 
same  as  from  Albert  and  Young  [2005].  Thus,  the  initial  flux 
is J{t  =  0)  =  cxp[-  (E  —  0.2)/0.1][sin(ao)  ~  sin(ao/.)],  where 
the  loss  eone  angle  ooi  =  5°  and  flux  j  is  related  to  phase 
space  density /by  j  ==flp‘.  Boundary  conditions  arc 

/L..,..  0.  (11) 


Ooo 


=  0. 


no 


n2) 


n 


E 


=  0. 


(13) 


where  Emm  “  0.2  MeV,  E^ax  =  5  MeV,  and  is  the 
momentum  corresponding  to  Emm  [Albert  and  Young, 
2005]. 

[16]  To  solve  the  equation  using  the  time-baekward  SDE 
method,  we  first  write  equation  (10)  in  the  form  of  (2): 


df  _ 

Ax(,(Jo  O-f  ,  M 


r  +  2 


Dn. 


Ot  p^  Ool  p  DaoOp  ^  <7/>- 


with 


h„„(l.  On, p)  = 


I  9  1  0  /CD„„A 

Gp  doo  \  p  )  G  Op  \  p  )' 


(15) 


(16) 
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“o  ao 

Figure  1.  D„o„o/p’.  and  |D,,()p|/p^:  inverse  time 

scales  in  units  of  s  ^  from  diffusion  coefficients  of  Albert 
and  Young  [2005].  Also  shown  is  the  sign  of  the  cross 
diffusion  eoeffieients.(Reprinted  from  Albert  and  Young 
[2005].) 

Thus,  the  two-dimensional  stochastic  process  defined  in 
equation  (5)  becomes 

dAois)  =  h^^(t-sMo>P)ds-haudW^  -hai2dfV2,  (18) 

dP(s)  =  bp(t-s,Ao,P)ds^  (721  dfVj  +(722  dW2.  (19) 

with  =  fVo  and  P(s=0)  =  p.  Then,  because  of  the 

Neumann  boundary  condition  at  ao  =  90°,  we  numerically 
reflect  Aq  with  respect  to  ao  =  90°  if  it  is  larger  than  90°. 
Here  components  of  the  matrix  a  are  defined  by 


(r\\(7\2 

fr\\  <721 

'^PiiapIP 

_<72I  <722  _ 

_<712  ^^22  _ 

[17]  In  this  work,  we  choose  <712  =  0  for  simplicity  and 
then  the  other  components  are 


■711  =  s/2D„„„Jp.  (21) 


£721  =  \/2Av,p/v/^  (22) 


£722  =  ^20^p  -  £7^1 ,  (23) 

where  we  have  used  the  fact  that  D^^oao  is  never  zero  in 
equation  (22). 

[18]  We  have  developed  a  2-D  SDE  code  to  solve  the 
diffusion  equation  (10)  where  SDEs  (18)  and  (19)  are 


integrated  using  the  Euler- Maruyama  method  [Kloeden 
and  Platen,  1992]: 

Aq  (>^0+1  )  —  •^o(‘^n)  "L  b(i^  [/  —  ,  /lo(5'rt),  /^(j'fl)] 

+  (S„)  Alfi  -f  (7i2(5„)  Alf2-  (24) 


^(^^«+l)  =  P(5n)  +/?/,[/ -5„,^o(.^«),/^(.V„)1  As 

-ha2i(s„)AfV,  -h  a22(s„)  AfV2.  (25) 

Here  AfV=  —  s„  N(0,  I),  where  iV(0,  I )  is  a  standard 
Gaussian  random  number  with  zero  mean  and  unit  variance 
generated  using  the  Box-Muller  algorithm  [Press  et  ai, 
2002].  Because  the  original  time-backward  SDE  method 
requires  fresh  samples  of  trajectories  for  every  different  (ao, 
p)  and  traces  trajectories  back  to  the  initial  condition  or  to  a 
boundary  every  time,  the  current  SDE  code  is  less  efficient 
when  solutions  on  many  grid  points  for  long  times  are 
needed.  Improving  the  efficiency  of  the  SDE  code  is  one  of 
tasks  in  our  future  work.  In  this  work,  we  mainly  want  to 
show  that  the  method  ean  be  used  to  solve  multidimensional 
diffusion  equations.  Results  from  the  SDE  code  are  compared 
with  those  of  Albert  and  Young  [2005]  in  section  3.2. 

3.2.  Comparisons  With  Results  of 
Albert  and  Young  [2005] 

[19]  Albert  and  Young  [2005]  solve  the  diffusion 
equation  (10)  by  first  transforming  to  new  coordinates 
which  diagonalize  the  diffusion  tensor  and  then  applying 
standard  finite  difference  methods  to  the  transformed  diffu¬ 
sion  equation.  The  bounee-averaged  diffusion  coefficients 
A»0o0»  ^o0/7»  and  Dpp  for  storm  time  chorus  waves  were 
calculated  at  Z.  =  4.5,  with  computational  methods  of  Albert 
[2005].  The  wave  model  used  to  calculate  diffusion  coef¬ 
ficients  is  described  by  Horne  et  al  [2005]  and  Albert  and 
Young  [2005];  the  wave  magnetic  field  is  given  by  ^5  = 
B^{u;)g^^{tan6),  where  the  wave  power  spectral  density 

and  the  wave  normal  angle  (tan^)  distribution  func¬ 
tion  ^^,(tan^)  are  truncated  Gaussian  functions  defined 
between  lower  and  upper  frequency  cutoffs  {u^ic  ^  ^  < 
u^uc)  and  wave  normal  angle  cutoffs  {Oi^  <  6  <  Oyc).  The 
latitudinal  distribution  of  the  waves  and  the  ratio  of  electron 
plasma  frequency  {fpe)  to  electron  cyclotron  frequency  (/^.^,) 
are  the  same  as  those  used  by  Home  et  al.  [2005]  and  Albert 
and  Young  [2005]  and  are  shown  in  Table  1 .  Similar  models 
were  used  by  Li  et  al.  [2007].  Up  io  n  =  ±5  resonance 
harmonies  were  included  in  the  calculation.  The  calculated 
diffusion  eoeffieients  AtOoo  are  proportional  to  (pAao)^  / 
A/,  as  from  Lyons  [1974a,  1974b]  and  are  divided  by  p^  to 
give  the  inverse  time  scales  plotted  in  Figure  I. 

[20]  Using  the  above  diffusion  coefficients  in  equation  (10), 
we  obtain  fluxes  for  £  =  0.5  MeV  and  2.0  MeV  electrons 
with  ao  ranging  from  6°  to  88°  with  1°  spacing  at  /  =  0.1 
and  I  day.  We  have  sampled  N  =  9000  trajectories  at  each 
ao  for  £  =  0.5  MeV  and  N  =  1 8000  trajectories  for  £  = 
2.0  MeV  with  dt  =  0.0004  day.  The  chosen  dt  gives  small 
relative  change  in  ao  and  £  per  step,  compared  with  scales 
of  the  diffusion  eoeffieients  and  initial  phase  space  density. 
Our  choices  of  N  and  dt  might  not  be  optimal,  and  choosing 
N  adaptively  is  probably  better  (G.  Cunningham,  personal 
communication,  2007).  Results  from  the  SDE  code  are 
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Table  1.  Latitudinal  Distribution  of  the  Waves  and  fpjfce  of  the 
Wave  Model  [Horne  et  ai,  2005]  Used  to  Calculate  Diffusion 
Coefficients 


Local  Time  Sector 

2300  0600  MLT 

0600-1200  MLT 

1200-1500  MLT 

fpjfce 

-^3.4  to  2.5 

~3.0  to  0.9 

-5.9  to  1.4 

Latitudinal  ranj^c 

0°  to  IS'* 

15°  to  35° 

10°  to  35° 

compared  with  those  of  Albert  and  Young  [2005].  Figure  2 
shows  the  comparisons  for  £  =  0.5  MeV  electrons  (Figure  2, 
top)  and  E  =  2.0  MeV  electrons  (Figure  2,  bottom),  with 
results  from  the  SDE  method  smoothed  using  a  six-point 
moving  window  average  in  ao  with  Aao  =  1^.  Within  small 
numerical  errors  associated  with  each  of  the  methods,  the 
two  sets  of  results  are  in  excellent  agreement,  and  they 
demonstrate  that  our  SDE  eode  is  able  to  suceessfiilly  solve 
the  bounee-averaged  pitch  angle  and  energy  diffusion 
equation. 

[2 1  ]  To  show  the  effeets  of  ignoring  off-diagonal  terms  on 
change  of  flux,  we  rerun  the  SDE  code,  setting  off-diagonal 
diffusion  coefficients  to  zero.  Results  are  shown  in  Figure  3 


Equotortol  pitch  angle  Qq  (degrees) 


Figure  2.  Comparisons  between  results  obtained  from  the 
SDE  method  (solid  lines)  and  the  Albert  and  Young  [2005] 
method  (dashed  lines)  for  (top)  E  =  0.5  MeV  and  (bottom) 
E  =  2.0  MeV  at  /  =  0.1  day  (blue  lines)  and  t  =  1  day  (red 
lines).  Here  black  lines  show  the  initial  condition. 


for  0.5  MeV  (Figure  3,  top)  and  2  MeV  (Figure  3,  bottom) 
eleetrons.  From  Figure  3  we  see  that  for  0.5  MeV  electrons, 
while  there  is  a  relatively  small  effeet  at  large  pitch  angles, 
ignoring  off-diagonal  terms  overestimates  electron  fluxes  at 
small  pitch  angles  by  a  factor  of  2  to  ^^5  at  /  =  1  day.  For 
2  MeV  eleetrons,  ignoring  off-diagonal  tenns  overestimates 
fluxes  by  a  factor  of  5  to  ~  1 0  at  /  =  1  day,  with  larger  errors 
at  smaller  pitch  angles.  Thus,  off-diagonal  terms  arc  more 
important  for  2  MeV  electrons.  Wc  emphasize  that  these 
results  are  for  the  Horne  et  al.  [2005]  wave  model,  and  wc 
note  that  the  peak  in  flux  of  2  MeV  electrons  near  30®  may 
be  related  to  the  cutoff  in  wave  power  at  35®  latitude  in  the 
Horne  et  al  [2005]  model  (see  discussion  in  section  4). 

3.3.  Effects  of  Parallel  Propagation  Approximation 

[22]  Summers  [2005]  and  Summers  et  al  [2007a,  2007b] 
have  derived  eyclotron  resonance  diffusion  coefficients  for 
field-aligned  waves,  where  only  Xhe  n  -  —1  resonanee  is 
ineluded  (heneeforth  denoted  by  D^*).  This  assumption  of 
parallel  propagation  greatly  improves  the  computation 
efficiency.  Bounce-averaged  D**  are  given  and  compared 
with  dif^sion  coefficients  obtained  from  the  PADIE  eode 


Equotoriol  pitch  ongle  ao  (degrees) 


Figure  3.  Fluxes  for  (top)  E  =  0.5  MeV  and  (bottom)  E  - 
2.0  MeV  at  /  =  0. 1  day  (blue  lines)  and  t  =  1  day  (red  lines) 
with  and  without  off-diagonal  diffusion  tenns.  Dashed  lines 
are  results  without  off-diagonal  diffusion  eoeffieients,  and 
solid  lines  are  results  with  off-diagonal  terms. 
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Figure  4.  Same  as  Figure  1,  except  that  diffusion 
coefficients  are  calculated  with  the  parallel  propagation 
approximation. 


[Glauert  and  Horne,  2005]  by  Shprits  et  al.  [2006].  In  the 
present  work,  we  also  calculate  D*  using  the  methods  of 
Albert  [2005]  with  the  same  wave  parameters  as  the  wave 
model  described  in  section  3.2,  except  that  ^lc  =  ^uc  ~  0. 
The  resulting  diffusion  coefficients  are  the  same  as  those 
obtained  from  the  PADIE  code  and  are  half  of  those  given 
by  Summers  et  al.  [2007a]  (this  factor  of  2  difference  is 
discussed  by  Albert  [2007]). 

[23]  Figure  4  shows  inverse  time  scales  from  diffusion 
coefficients  with  the  parallel  wave  approximation.  Com¬ 
pared  with  Figure  1 ,  wc  see  that  the  general  behavior  of  d'' 
is  quite  good,  with  larger  differences  for  £  >  1  MeV 
electrons.  The  off-diagonal  terms  of  d"  arc  worse 
approximations  than  the  diagonal  terms,  with  details 
discussed  by  Albert  [2007]. 

[24]  To  compare  effects  of  D*'  with  fiilly  calculated 
diffusion  coefficients  D,  we  solve  equation  (10)  for  0.5  MeV 
and  2  MeV  electrons  using  the  following  four  sets  of 
diffusion  coefficients:  (1)  D**,  (2)  diagonal  tenns  of  d" 
(hereinafter  referred  to  as  \  (3)  D,  and  (4)  diagonal  terms 
of  D  (hereinafter  referred  to  as  Dj).  Results  arc  shown  in 
Figures  5  7. 

[25]  Figure  5  (top)  shows  the  comparison  between  fluxes 
calculated  using  and  Dj  for  0.5  MeV  electrons.  We  see 
that  results  from  Dj  agree  very  well  with  with  slight 
differences  for  Qq  greater  than  about  40°.  Figure  5  (bottom) 
shows  the  same  comparison  for  2.0  MeV  electrons  from 
which  we  see  that  the  flux  from  dJ  is  smaller  than  that  from 
Dj  by  up  to  ~5  orders  of  magnitude  at  low  ao  (<1 5°)  at  t  = 
1  day.  This  behavior  occurs  because  Djj  underestimates 
energy  diffusion  coefficients  for  high-energy  particles  at 
small  pitch  angles,  where  n  ^  —\  resonances  also  make  a 
significant  contribution.  Thus,  d\  produces  larger  differ¬ 


ences  in  fluxes  for  2  MeV  electrons  than  0.5  MeV  at  small 
ao  compared  with  Dj. 

[26]  Figure  6  shows  comparisons  between  fluxes  calcu¬ 
lated  using  and  D  for  0.5  MeV  electrons  (Figure  6,  top) 
and  2  MeV  electrons  (Figure  6,  bottom).  Figure  6  (top) 
shows  that  dJ/  overestimates  increase  of  flux  at  small  pitch 
angles  for  0.5  MeV  electrons,  which  is  expected,  because  d|1 
yields  very  similar  flux  increases  as  Dj  for  0.5  MeV 
electrons.  For  2.0  MeV electrons,  fluxes  from  D^are  smaller 
than  that  from  D  for  ao  ^  18°  and  larger  for  ao  ^  1 8°  at  /  = 
1  day  (where  the  difference  can  be  about  1  2  orders  of 
magnitude). 

[27]  Fluxes  calculated  from  D"  and  D  (i.e.,  with  off- 
diagonal  terms  included)  for  0.5  MeV  and  2  MeV  electrons 
are  shown  in  Figure  7  (top)  and  Figure  7  (bottom), 
respectively.  Reasonable  agreement  between  D  and  D 
fluxes  is  obtained  for  ao  ^  50°,  but  significant  differences 
occur  at  smaller  pitch  angles.  For  0.5  MeV  electrons,  d" 
underestimates  increases  of  flux  at  /  =  1  day  by  approxi¬ 
mately  an  order  of  magnitude  for  ao  <  20°.  For  2.0  MeV 
electrons,  behavior  of  D*'  is  worse  at  /  =  1  day.  We  see  from 
Figure  7  (bottom)  that  D''  underestimates  increases  of  flux 


Equotorlol  pitch  angle  Oq  (degrees) 


Figure  5.  Comparisons  between  results  obtained  from 
dififtision  coefficients  (dashed  lines)  and  Dj  (solid  lines) 
for  (top)  E  =  0.5  MeV  and  (bottom)  E  =  2.0  MeV  at  /  =  0.1 
day  (blue  lines)  and  /  =  1  day  (red  lines). 
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Equotoriol  pitch  ongle  a©  (degrees) 


Figure  6.  Comparisons  between  results  obtained  from 
diffusion  coefficients  Dj  (dashed  lines)  and  D  (solid  lines) 
for  (top)  E  =  0.5  MeV  and  (bottom)  E  -  2.0  MeV  at  /  = 
0.1  day  (blue  lines)  and  t  =  1  day  (red  lines). 

by  ~  1  4  orders  of  magnitude  for  10®  :S  ao  ^  35®.  Thus, 
the  approximation  of  parallel  propagation  produees  larger 
differenees  at  small  piteh  angles  for  higher-energy  partieles, 
especially  when  off-diagonal  terms  are  included. 

4.  Summary  and  Discussion 

[2s]  In  this  work  a  new  eode,  based  on  the  mathematieal 
theory  of  expressing  solutions  of  diffusion  equations  in 
terms  of  related  stochastic  processes,  has  been  developed 
for  solving  multidimensional  radiation  belt  diffusion 
equations.  Two  examples  are  used  to  show  its  applieations. 

[29]  First,  we  apply  the  SDE  eode  to  a  bounee-averaged 
piteh  angle  and  energy  diffusion  equation  and  obtain 
exeellcnt  agreement  with  a  previously  developed  method 
[Albert  and  Young,  2005].  We  also  eonfirm  that  ignoring 
off-diagonal  terms  in  the  diffusion  equation  overestimates 
inerease  of  flux,  espeeially  at  small  piteh  angles,  at  /  =  1  day 
(by  a  factor  of  2  to  ^^5  for  0.5  MeV,  and  5  to  ~  1 0  for  2  MeV 
electrons)  using  the  Albert  and  Young  [2005]  diffusion 
coefficients. 

[30]  Second,  by  solving  the  bounce-averaged  pitch  angle 
and  energy  diffusion  equation  using  fully  calculated  diffu¬ 


sion  coefficients  D  [Albert  and  Young,  2005]  and  coeffi¬ 
cients  with  the  parallel  propagation  approximation  D** 
[Summers,  2005;  Summers  et  ai,  2007a,  2007b],  both 
calculated  using  the  chorus  wave  model  of  Horne  et  al. 
[2005],  we  show  that  diagonal  difflision  coefficients  of  D** 
agree  well  with  those  of  D  only  for  low-energy  particles 
(e.g.,  £'=0.5  MeV).  For  high-energy  electrons,  the  difference 
between  the  diagonal  terms  of  D  '  and  D  produces  large 
differences  in  fluxes  at  some  pitch  angles  (difference  of  up 
to  5  orders  of  magnitude  for  2  MeV  electrons  at  Oo  15®,  at 
t  =  1  day).  By  including  off-diagonal  diffusion  coefficients 
in  our  calculation,  we  show  that  the  off-diagonal  terms  of 
can  produce  differences  in  fluxes  of  4  orders  of 
magnitude  for  2  MeV  electrons  at  /=!  day.  A  discussion 
of  the  details  of  different  diffusion  coefficients  and  another 
approximation  for  a  full  calculation  of  diffusion  coefficients 
are  presented  by  Albert  [2007]. 

[31]  Note  that  the  above  conclusions  on  the  magnitude 
and  location  of  differences  that  occur  by  omitting  off- 
diagonal  terms  and  assuming  parallel  propagating  waves 
are  very  likely  to  be  dependent  on  the  wave  model  used.  For 
example,  a  different  latitudinal  distribution  of  wave  power 
may  result  in  different  diffusion  coefficients  and  thus 


Equotoriol  pitch  angle  (degrees) 


Figure  7.  Comparisons  between  results  obtained  from 
diffusion  coefficients  D''  (dashed  lines)  and  D  (solid  lines) 
for  (top)  £  =  0.5  MeV  and  (bottom)  £  =  2.0  MeV  at  t  - 
0.1  day  (blue  lines)  and  t  =1  day  (red  lines). 
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Figure  Al.  Local  effects  of  ignoring  off-diagonal  terms. 
Lines  are  contours  of  particle  numbers.  Particles  are 
released  from  ao  =  30°  for  £  =  1 .0  MeV  and  ao  =  50°  for 
E  =  3  McV,  (left)  Off-diagonal  terms  are  kept.  is 

positive  at  aQ  =  30°,  £  =  1.0  MeV  and  negative  at  qq  =  50°, 
£  =  3  MeV,  (right)  Off-diagonal  terms  are  set  to  zero. 


different  conclusions.  The  sensitivity  of  our  results  to  wave 
models  needs  further  study.  However,  before  such  work  is 
done,  it  is  safer  to  include  both  off-diagonal  terms  and 
oblique  waves  in  calculations  of  electron  flux, 

[32]  The  SDE  method  is  less  efficient  when  solutions  on 
many  grid  points  are  desired.  However,  when  parallel 
computers  are  available,  computation  time  can  be  greatly 
reduced  because  of  high  parallelization  effleieney.  General¬ 
ization  to  3-D  including  pitch  angle,  energy,  and  radial 
diffusion  is  straightforward.  The  SDE  method  is  very 
promising  for  providing  new  insights  into  the  relative  roles 
of  local  acceleration  and  radial  diffusion  as  acceleration 
mechanisms  and  the  importance  of  pitch  angle  diffusion  as  a 
loss  process. 


Appendix  A:  Time-Forward  SDE  Method 

[33]  To  use  the  time-forward  SDE  method,  we  first  set 
F  =  Gf  and  write  the  bounee-averaged  pitch  angle  and 
energy  diffusion  equation  (10)  in  the  following  form: 


Oal 


+  2 


OaoOp 


(Al) 


where  and  bp  are  defined  in  equations  (16)  and  (17). 
Thus,  the  time-forward  stochastic  differential  equations 
corresponding  to  equation  (Al)  are  [  Alanko-Huotah  et  ai, 
2007;  Yamada  et  al,  1998;  Qin  et  al.,  2005] 


dA^it)  =  b^,(t,A^,P)dt^ay^dW^  ^a^2dW2.  (A2) 


dP{t)  =bp{t,AQ.P)dt  +  (j2\dW^  -\-(J22dW2.  (A3) 

where  components  of  the  matrix  cr  are  also  defined  by 
equations  (21)  (23), 


[34]  Equations  (A2)  and  (A3)  are  solved  to  give  changes 
of  particle  coordinates  (ao,  p).  Thus,  after  a  given  time 
period,  the  distribution  of  electrons  can  be  obtained.  Here 
we  choose  a  time  period  short  enough  to  ignore  boundary 
effects.  To  explore  local  effects  of  off-diagonal  diffusion 
coefficients  on  distributions  of  particles,  we  release  9000 
particles  from  ao  =  30°,  £  =  1  MeV,  where  D^^^p  is  positive, 
and  ao  =  50°,  £=  3  MeV,  where  D^^o/?  is  negative.  We  obtain 
the  distribution  of  particles  shown  in  Figure  A 1  after  f  =  0.06 
day  for  £  =  3  MeV  and  /  =  0.0 1  day  for  £  =  1  McV.  Wc  also 
turned  off-diagonal  diffusion  coefficients  on  and  off  to  show 
local  effects  of  ignoring  off-diagonal  terms.  Figure  Al  (left) 
has  D(>o/7  ^  0,  and  Figure  Al  (right)  has  AvO/?  ^  0.  We  see 
from  Figure  Al  that  without  » the  local  distribution  of 
particles  has  a  shape  of  an  ellipse,  while  with  D,»o/7  this 
ellipse  is  tilted,  and  the  tilt  direction  is  determined  by  the 
sign  of  D,.,i)p.  With  positive  (as  for  the  ao  =  30°,  £  =  1 
MeV  ease)  the  ellipse  tilts  clockwise,  and  with  negative 
(ao  =  50°,  £  =  3  MeV),  the  ellipse  tilts  counterclockwise. 
These  results  are  consistent  with  previous  analytical  results 
using  Green  functions  [Albert  and  Young,  2005]. 
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